function [samples, labels] = gmmrnd(priors, mmu, ssigma, NUM_SAMPLES, PERMUTE_SAMPLES)
% generate labeled testdata 
% [data, labels] = genTestData(priors, mmu, ssigma)
%
% priors ... length = size(mmu,1) or size(mmu,1)+1 for uniformly distributed noise
% mmu ... (NUM_KERNELSxDIM)
% ssigma ... (DIMxDIMxNUM_KERNELS) or 
%            (NUM_KERNELSxDIM²) for cholesky factorizations
% PERM_SAMPLES ... (optional) if true, samples are permuted
    
    
    if nargin < 5, PERMUTE_SAMPLES = false; end
    
    if sum(priors)-1 > eps
        fprintf('\n\n *EE* sum of priors <> 1, aborting. *EE*\n');
        return;
    end
    
    [NUM_KERNELS DIM] = size(mmu);
    
    num = round(NUM_SAMPLES*priors);
    diff = NUM_SAMPLES - sum(num); % make sure to gen NUM_SAMPLES samples
    ind_k = ceil(NUM_KERNELS*rand(abs(diff),1));
    num(ind_k) = num(ind_k) + sign(diff);
    
    disp(sprintf('\n *II* generating %d %d-dimensional samples *II*\n', ...
                 sum(num), DIM));

    samples = zeros(NUM_SAMPLES,DIM);
    labels = zeros(NUM_SAMPLES,1);
    
    ind1 = 1;
    for k = 1:NUM_KERNELS
        ind2 = ind1+num(k);
        R = reshape(ssigma(k,:),DIM,DIM);
        samples(ind1:ind2-1,:) = mvnrnd(mmu(k,:),R'*R,num(k));
        labels(ind1:ind2-1) = k;
        ind1 = ind2;
    end

    if PERMUTE_SAMPLES
        disp(' permuting samples ...');
        P = spalloc(NUM_SAMPLES, NUM_SAMPLES, NUM_SAMPLES);
        perm_ind = NUM_SAMPLES * [0:NUM_SAMPLES-1] + randperm(NUM_SAMPLES);
        P(perm_ind) = 1;
        samples = P*samples;
        labels = P*labels;
        clear P perm_ind;
    end

    %gmmplot(samples, mmu, ssigma);
    
end
